
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!c
!c    routine to invert a matrix by Gaussian elimination
!c     Ainv=inverse(A)    A(n,n)   Ainv(n,n)
!c     Note: sizes maxed at 64 x 64 - re dimension for larger matrices;
!c           D matrix is extended.
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       subroutine invert( A, Ainv, n )
       implicit double precision (a-h,o-z)
       double precision A(n,n),Ainv(n,n)
       double precision D(n,2*n)


!c  initialize the reduction matrix
       n2 = 2*n
       do  i = 1,n
        do  j = 1,n
         D(i,j) = A(i,j)
         d(i,n+j) = 0.
        enddo
        D(i,n+i) = 1.
       enddo

!c  do the reduction 
       do 3 i = 1,n
        alpha = D(i,i)
        if(alpha .eq. 0.) go to 300
        do 4 j = 1,n2
         D(i,j) = D(i,j)/alpha
 4      continue
        do 5 k = 1,n
         if((k-i).eq.0) go to 5
         beta = D(k,i)
         do 6 j = 1,n2
          D(k,j) = D(k,j) - beta*D(i,j)
 6       continue
 5      continue
 3     continue

!c  copy result into output matrix
       do 7 i = 1,n
        do 8 j = 1,n
         Ainv(i,j) = D(i,j+n)
 8      continue
 7     continue


       return
 300   print *,'*** ERROR: Singular matrix *** A = ', A
       Ainv=0.0d0
       return
       end
